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Abstract 

The slope of the coexistence line of the liquid-liquid phase transition (LLPT) can be positive, 
negative, or zero. All three possibilities have been found in Monte-Carlo simulations of a mod- 
ified spherically symmetric two-scale Jagla model. Since the liquid-liquid critical point (LLCP) 
frequently lies in a region of the phase diagram that is difficult to access experimentally, it is 
of great interest to study critical phenomena in the supercritical region. We therefore study the 
properties of the Widom line, which is defined in the one-phase region above the critical point as 
the locus of maximum correlation length as function of the ordering field at constant thermal field. 
Asymptotically close to the critical point, the Widom line coincides with the loci of the response 
function extrema, because all response functions can be asymptotically expressed as functions of 
the diverging correlation length. We find that the method of identifying the Widom line as the 
loci of heat capacity maxima becomes unfruitful when the slope of the coexistence line approaches 
zero in the T-P plane. In this case the specific heat displays no maximum in the one-phase region 
because for a horizontal phase coexistence line, according to the Clapeyron equation, the enthalpy 
difference between the coexisting phases is zero, and thus there can be no contribution to enthalpy 
fluctuations from the critical fluctuations. The extension of the coexistence line beyond the critical 
point into the one-phase region must in this case be performed using density fluctuations; the 
line of compressibility maxima remains well defined, though it bifurcates into a symmetrical pair 
of lines. These findings agree well with the linear scaling theory of the LLCP by Anisimov and 
collaborators. 
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I. INTRODUCTION 



The liquid-liquid phase transition (LLPT), defined as a transition between two liquid 
states of different densities, called low density liquid (LDL) and high density liquid (HDL), 
has received considerable attention not only due to its rarity in nature, but also due to its 
importance in our fundamental understanding of the liq uid state of matter Q-Q. The LLPT 
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ionic 12j, molecular 



13J, and 



was observed in many systems, such as elemental 
covalent 14j liquids. In some cases, the LLPT can terminate at a liquid-liquid critical point 



(LLCP). Systems such as liquid water, silicon, silica, and germanium, possess analogous 



thermodynamic and dynamic anomalies |3j, [10J, [15-20]. The critical phenomena near the 
LLCP are of crucial importance for the understanding of the anomalous properties in these 

0,y,r 



21] . However the detection of the LLCP or the LLPT can be difficult 



systems 



due to the fact that in many cases, the LLCP is deeply buried in the supercooled region, 
where crystallization may occur before we reach the LLCP Q|. For example, in the case 
of water, it has been hypothesized that the putative LLCP is the cause of water anomalies 
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231 ]. but the existence of a LLCP for bulk water in the deeply supercooled region has 



not been directly verified by experiments due to crystallization, even though indications of 
the existence of the LLCP have been found both in pressure induced melting experiments 



241 ] and in nanoconfined water 201 ] . 



According to scaling theory, asymptotically near the critical point all response functions 
can be expressed in terms of the correlation length 25J] . Different response functions diverge 
at the critical point, and display maxima in the one-phase region along constant pressure P 



paths or constant temperature T paths 
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261 ] . The loci of different response function 



maxima in the T-P plane are different, but they converge in the vicinity of the critical point 
to a single line, called the Widom line (3), [21], 27 1. Theoretically, Widom line is defined as 
the locus of maximum correlation length £ as function of ordering field at constant thermal 
field in the one-phase region. Approaching the critical point, the magnitude of the response 
functions increases, and becomes infinite at the critical point. This fact provide a new way 
of locating the critical point: instead of locating the critical point through the coexistence 
line below the critical temperature T c , we may locate the critical point through the Widom 
line in the one-phase region from the higher temperature side by studying the behavior of 
the loci of response function maxima [3]. Thus it is important to find a general model system 
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with an accessible LLCP which would permit detailed examination of the response functions 
in the vicinity of the LLCP. 

The Jagla model of liquids is a simplified model consisting of particles interacting via a 



spherically symmetric two-scale potential with both repulsive and attractive ramps 



28 
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29] . With a special choice of parameters, the Jagla model has an accessible LLCP above 



the melting line {3], allowing us to explore the behavior near the LLCP in equilibrium liquid 
states. In this case, the coexistence line between LDL and HDL is positively sloped, which 
means that when cooled down along the same isobar, the system changes from LDL to HDL. 



This behavior is opposite to that of water, where experiments |20| and simulations 21] show 
that the coexistence line might be negatively sloped, and an isobaric cooling path transforms 
the system from HDL to LDL jj. 

Gibson and Wilding found that by changing the parameters of the Jagla potential it is 
possible to reduce the slope of the coexistence line to zero [30] . In this paper we use modified 
Jagla models to investigate the behavior of the Widom line as the slope of the coexistence 
line changes from positive to horizontal. In Sec. II we introduce the modified Jagla model 
and the simulation method. In Sec. Ill we present our simulation results. In Sec. IV we 
compare our simulation results to the linear scaling theory of the critical point. In Sec. V 
we further investigate the relationship between the LLCP, the Widom line, and the glass 
transition for systems with different coexistence line slopes. In Sec. IV we summarize our 
study. 



II. MODEL 



Here we study the two length-scale Jagla model with both repulsive and attractive ramps 
. In this model, particles interact with a spherically symmetric pair potential 



U(r) 



(1) 



oo r < a 

U A + (U R - U A ){b - r)/(6 - a) a < r < b 

U A (c — r) / (c — b) b<r<c 
r > c 

where a is the hardcore distance, b is the soft-core distance, and c is the long-distance cutoff 
[Fig. [I]; Ua = —Uq is the minimal potential energy reached at soft-core distance r = b, and 
Ur is the potential energy at the top of the repulsive ramp at hardcore distance r = a. 



We implement a family of Jagla potentials with different parameters, simultaneously de 



creasing b and c — essentially following the Gibson- Wilding procedure [30j, the only difference 
being that we keep U a constant. The parameters of different models are presented in Table[H 
We perform discrete molecular dynamics (DMD) simulations by discretizing the ramp 
into a series of step functions. The discrete Jagla potentials are 

oo r < a 

Ur a < r < a + |An 

Un-kAUt a + (k - |)An < r < a + (k + 5) An, 1 < k < n x - 1 

U A b- |An < r < b + §Ar 2 

U A + kAU 2 b+{k-\)Ar 2 <r < 6 + (A; + |)Ar 2 , 1 < fc < n 2 - 1 
r > c - |Ar 2 

where ri\ = 60 and n 2 = 20, An = (b — a)/ni, AUi = (Ur — Ua)/tii, and Ar 2 = (c — b)/n 2 , 
AU 2 = U /n 2 . 

We use a as the unit of length, particle mass m as the unit of mass, and Uq as the 



U k (r) 



(2) 



unit of energy. The simulation time t is measured in units of a^m/Uo, temperature T in 
units of Uo/ks, pressure P in units of Uo/a 3 , density p = N/L 3 in units of a -3 , isothermal 
compressibility in units of a 3 /Uq, and isobaric specific heat Cp in units of ks- 

Our results are based on simulations of a liquid system of N = 1728 molecules with 
periodic boundary conditions. Constant volume-temperature (NVT) and constant pressure- 
temperature (NPT) simulations are implemented in this study. 

The temperature of the system is controlled by rescaling the velocities of all particles in 
the NVT simulations so that the average kinetic energy per particle approaches the desired 
value 3KrT /2, where T is the temperature of the thermostat, 

T' = T(l- K T T t )+T Q K T T U (3) 

where kt = 0.2 [y/m/Uo/a] is a heat exchange coefficient, Tj is the time interval between two 
successive rescaling, T' is the new temperature, and T is the average temperature during 
the time interval r t . We select r t as the time during which N collisions between particles 
occur. 

For the NPT simulations, the Berendsen barostat algorithm rescales the coordinates fj 
and box dimensions L after each r v unit of time, 

r'j = r 3 + r 3 K p (P - P ) (4) 



L = L + Lk p {P-P ), (5) 

where Pq is the desired pressure, P is the average pressure during time interval r v = 1000r t , 
and k p = 0.02[a 3 /£/o] is the rescaling coefficient. 

We perform cooling or heating simulations at a constant cooling/heating rate, q = 
AT /At, where the T decreases/increases by AT over time At. We measure q in units 
of g = x/U$/(ma 2 k 2 B ). 



III. SIMULATION RESULTS 



A. Coexistence line 



We first explore the phase diagram of each model with different b/a via slow cooling using 
constant volume simulations. The LLCP corresponds to the highest temperature crossing of 
isochores in the T-P phase diagram. The temperature of maximum density (TMD) line is the 
locus of state points at which pressure reaches a minimum along each isochore as a function 
of T {27]]. Figure EJ^a) shows that the LLCP monotonically shifts to lower temperature and 
higher pressure as b/a decreases. The region of the density anomaly (the region bounded 
by the TMD line) expands and shifts together with the LLCP to higher pressures as b/a 
decreases. This behavior coincides with the behavior reported in Ref. 30[. (The numerical 
differences in P and T arise from the fact that we define the unit of energy as Uq = —Ua, 
while Ref. [30] uses Uo = Ur/0.69.) When b/a < 1.59 the system crystallizes spontaneously 
near the LLCP within a short simulation time, and we are not able to obtain the LLCP and 
coexistence line. 

We calculate the slope of the LLPT coexistence line dP/dT for systems with different 
b/a [Fig. Et^b)] by the Maxwell construction. One can see that the slope decreases from 
positive to approximately zero as b/a decreases to 1.59, in agreement with Ref. 30[. For 
large b/a, the LLCP lies clearly above the density anomaly region bounded by the TMD 
line, corresponding to the case of dP/dT > 0, while for b/a = 1.59, the LLCP lies on the 
TMD line corresponding to the case of dP/dT = [Fig. 12fa)|. Theoretically, if dP/dT < 0, 



the LLCP should be inside the density anomaly region 
scaling theory, which we will address in detail in Sec. VI. 



31], as is confirmed by the linear 



6 



B. Widom line 

For the first order phase transition, the order parameter, entropy, or density is discon- 
tinuous on crossing the coexistence line. At the critical point, where the coexistence line 
terminates, the critical fluctuations of Cp and diverge, and show maxima in the one 
phase region close to the critical point. In this section, we study the behavior of Cp maxima 
and Kt maxima lines near the LLCP for models with different coexistence line slopes. 

1. Isobaric specific heat Cp 

We first explore the behavior of the Cp maxima line for modified Jagla models with 
different b/a. For b/a > 1.62, the coexistence line has a positive slope [Fig. [2(b)] . Upon 
cooling along constant pressure above the critical pressure P c in the one-phase region, we 
observe peaks in Cp for the cases b/a = 1.72, 1.70, 1.68, and 1.65 [Fig. [3(a-d)]. As the LLCP 
is approached the magnitude of the Cp peaks increases, and at the LLCP they diverge in 
an infinite system. When P < P c , we observe a continuous increase in Cp without any peak 
when the coexistence line is crossed. Thus we can locate the LLCP by locating the terminal 
point of the Cp maxima line. 

We note that the Cp peaks move toward higher T at higher P, indicating a positively 
sloped line [Fig. [3(a,b)]. This is consistent with the fact that the Widom line is the ex- 
tension of the coexistence line into the one-phase region, and for these values of b/a the 
coexistence line is positively sloped. However the slopes of the Cp maxima lines increase as 
the slopes of the coexistence lines decrease and eventually, at b/a = 1.65, the Cp maxima 
line becomes nearly vertical, clearly showing that the Cp maximum is no longer serving its 
original purpose, as will be explained in Sec. IV below. For b/a = 1.62, Cp monotonically 
increases without showing any peak for pressures P > P c , except at the highest pressure 
studied P = 0.800 [Fig. [3(e)]. 

When b/a = 1.59 with a horizontal coexistence line, Cp monotonically increases with 
decreasing T along a constant pressure path both below and above P c [Fig. [3(f)]. When 
b/a = 1.62, there are no Cp maxima in the equilibrium region with T > T c , but Cp behaves 
symmetrically either below or above P c . 

We plot the lines of equal Cp for two extreme cases, b/a = 1.72 with a positively sloped 
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coexistence line [Fig. Ufa)] and b/a = 1.59 with a horizontal coexistence line [Fig. Hfb)]. 
When b/a = 1.72, the lines of equal Cp form loops in the T-P plane and cross the Cp 
maxima line at their maxima points. The locus of Cp maxima extends the coexistence line 
into the one-phase region in the vicinity of the critical point. Then it sharply turns upwards 
to higher pressures and becomes approximately vertical. For b/a = 1.59, there are no Cp 
maxima. The equal Cp lines extend away from the critical point symmetrically without any 
loops. At low T, we reach the simulation limit due to either crystallization for P < P c or 
due to entering a glassy state for P > P c , where no equilibrium results can be obtained for 
the analysis. 

We note that the magnitude of Cp drops significantly when the coexistence line is hori- 
zontal with b/a = 1.59, compared to when b/a = 1.72. This is because, when the coexistence 
line is horizontal, the difference in enthalpy H between LDL and HDL is zero according to 
the Clapeyron equation of thermodynamics, 

— = — (6) 
dT TAV [ ' 

In this case, the enthalpy fluctuations that determine the magnitude of the specific heat 
gain no strength from the critical fluctuations, except from the weak PV term. 



2. Isothermal compressibility Kt 

Figure |5] shows the behavior of Kt above and below P c and contrasts it with that of Cp. 
For b/a > 1.62, when the coexistence line slope is positive, Kt shows maxima both above 
and below P c . For P > P c in the one-phase region, similar to Cp, the Kt peaks become more 
prominent as the LLCP is approached [Fig. |5^a-e)]. For a finite system, Kt diverges at the 
LLCP. For P < P c , we observe a second set of Kt peaks with much lower magnitudes. For 
b/a = 1.59 with a horizontal coexistence line, Kt behaves symmetrically above and below 
P c , with equal magnitudes of the maxima [Fig. [5(f)]. 

Figure [6] shows the loci of Kt maxima for both b/a = 1.72 with a positive coexistence 
line slope and b/a = 1.59 with a horizontal coexistence line. For b/a = 1.72, the values of 
the Kt maxima at P > P c , which corresponds to the critical fluctuations and originates 
from the LLCP, have a much larger magnitude than the values of Kt maxima at P < P c . 
The second Kt maxima line at P < P c corresponds to the approach to the LDL spinodal 
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and terminates at the lowest pressure point of the LDL spinodal where the TMD line also 



terminates 



22 



32] . Furthermore, the P < P c Kt maxima line also crosses the TMD line at 



the point of its maximal temperature 31] . Similar to Cp, the lines of equal Kp form loops, 
and cross at their maximal pressure points with the P > P c and P < P c Kt maxima lines. 
For b/a = 1.59 [Fig. [5], the lines of equal Kt form symmetric (with respect to P c ) loops 
around the critical point. Both Kt maxima lines merge at the LLCP. 

In the case of b/a = 1.72, from both Fig. H] and Fig. [6], we can identify the Widom line 
as the overlapping segment of the Cp maxima and the Kp maxima lines, which extends 
the coexistence line into the one-phase region in the vicinity of the LLCP. In contrast, for 
b/a = 1.59 with a horizontal coexistence line, the Cp maxima line disappears, where the 
specific heat Cp can no longer be a good representative for critical fluctuations. Indeed, there 
is no enthalpy difference between the two coexisting phases, so there can be no contribution 
to enthalpy fluctuations from the critical fluctuations. However, the density fluctuations 
remain well defined, with two Kt maxima lines above and below the critical pressure, both 
associated with the critical fluctuations. In the vicinity of the critical point, the two Kt 
maxima lines merge together and can be used to locate the critical point from measurements 
obtained in the supercritical region only. 



IV. COMPARISON WITH THE LINEAR SCALING THEORY OF THE LIQUID- 
LIQUID CRITICAL POINT 

To explain the change of behavior of the lines of response function extrema for the hor- 
izontal coexistence line case, we adapt the linear scaling theory of the liquid-liquid critical 
point developed by M. A. Anisimov and collaborators [33] . 

The field-dependent thermodynamic potential ip can be considered a universal function of 
two scaling fields: "ordering" h\ and "thermal" \\2- Near the critical point ip can be written 

where the critical exponents have the values for the Ising universality class, a = 0.109, 
P = 0.326, and 7 = 1.239 J33 ] . 



Since our focus is on the immediate vicinity of the liquid-liquid critical point, we neglect 
the curvature of the coexistence line and the background contribution to the response func- 
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tions. We assume the scaling fields are linear analytical combinations of physical fields, the 
pressure P, and the temperature T, 

h = a x AP + a 2 AT, (8) 

h 2 = hAT + b 2 AP. (9) 

with AP = (P - P c )/(p c RT c ) and AT = (T — T c )/T c , where the subscript "c" here and 
below indicates the critical parameters, and aj and 6, are system-dependent coefficients. 

We introduce a tuning parameter into the theory, and we use it to change the slope of 
the coexistence line by defining 

a\ = cos 0, a 2 = — sin 0, 

b\ = cos0, 6 2 = sin0. (10) 

Then the slope of the coexistence line is 

dP 

V =*>»*■ (ID 

The critical (fluctuation-induced) parts of the dimensionless isobaric specific heat and 
isothermal compressibility are expressed through the scaling susceptibilities j33| 

<cw = n% )p 

= T(sin 2 0xi - sin 20xi2 + cos 2 0X2), (12) 

= -i( d -t )f 

V dP 

= ^- [cos 2 0xi - sin 20xi2 + sin 2 0x 2 ] , (13) 

(ctp) = — ( — -)p 
V dT 

= sin20xi - cos 20xi2 - ^ sin20x 2 )- (14) 

where f = T/T c , P = P/p c RT c . 

The scaling fields and scaling susceptibilities can be written as functions of the "po- 
lar" variables r and 8, and two constants a and k, which can be obtained by fitting the 
experimental data, 

h x =ar p+1 6{l -6 2 ), (15) 
10 



h 2 = r(l-ti l l ), (16) 

Xi = -r-'c^e), (17) 
a 

X12 = kr?- l c 12 (e), (18) 

X2 = akr- a c 2 (6), (19) 



where 



c 1 (9) = [l-b 2 9 2 (l-2(3)]/c (9), (20) 

cu(9) = (39[l - 5 - 9 2 (3 - S)]/c (9), (21) 

C2(6) = [(1 - a)(l - 39 2 )s(9) - 2s 2 f359 2 (l - 9 2 )]/c (9), (22) 

co(0) = (1 - 36» 2 )(1 - b 2 9 2 ) + 2/35b 2 9 2 (l - 9 2 ). (23) 

Here b = ^(7 - 2^)/7(l - 2/3) = 1.16679, 6 = j/P+l = 4.80061, and s(0) are known 
functions, 

s (0) = So + S2 2 , (24) 
So = -(2 - a)/ , (25) 
S2 = -(2-a)6 2 (l-2/3)/ -7/ 2 , (26) 

with 

g(g _ 3) _ &2q7 
70 26 4 (2-a)(l-a)a' 1 J 

_ 3) -^(1-2/3) 
/2 " 26>(1 - a)a " (28) 

Since all response functions are directly propositional to k, its actual value is irrelevant 
for our study. The value of a determines the strength of the ordering field. We found that for 
large values of a, the overlapping segments of Cp and Kt maxima lines are shorter, which 
more adequately models the behavior of our simulations. 

We study the behavior of Cp and Kt in systems with different coexistence line slopes. 
Figure [7] shows a plot of the loci of Cp and Kt maxima in the T-P plane. Note that both 
the locus of Cp maximum and the Kt maxima line with the higher magnitude originate 
from the LLCP, and coincide with each other close to the LLCP before they separate. The 
locus of Cp maximum bends towards a lower temperature than that of the Kt maximum 
line. As the slope of the coexistence line decreases, the overlapping segments of the loci of 
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Cp and Kt maxima shorten. When the coexistence line is horizontal, the loci of Cp and 
Kt maxima separate. 

We then examine the systems with the slope of coexistence line very close to zero with 
small <p as shown in Fig. |HJ We see that the loci of Kt maxima are similar, with two Kt 
maxima lines approaching the LLCP from T > T c , but the Cp maxima line deviates from 
the Kt maxima lines and its slope becomes more vertical, as the slope of the coexistence 
line approaches zero. For the horizontal coexistence line, a second Cp maxima line emerges, 
and converges with the first Cp maxima line at T < T c . Both of the Cp maxima lines enter 
the critical point almost horizontally from T < T c , while the Kt maxima lines enter from 
T > T c horizontally. Thus in this case, close to the LLCP at T > T c , no Cp maxima can 
be found. According to the linear scaling theory, the actual Cp maxima may exist below 
T c where, in our simulations, they are buried in the crystallization region or in glassy states 
with no equilibrium data for analysis. Further, there is no convergent behavior of the Cp 
maxima and Kt maxima near the LLCP in the horizontal coexistence line case. The Widom 
line, approximated by the locus of density fluctuation maxima, is the convergent loci of the 
Kt maxima, while the identification of the Widom line using the Cp maxima is no longer a 
fruitful method. 

We also find out that when the slope of the coexistence line is positive, the contours 
of Cp and Kt go around the LLCP, with the pressure maxima following the extension of 
the coexistence line into the one-phase region [Fig. [9^a-f)]. When the slope of coexistence 
line is horizontal, there are symmetric contours above and below the critical pressure. This 
matches well with what we find in our simulation [Figs. H] and [6] . 

V. GLASS TRANSITION 

By decreasing the b/a ratio, the LLCP is pushed into a metastable region with respect to 
crystallization, where the system is close to the glass transition (GT). We hence investigate 
the relationship between LLPT and GT in systems with different slopes of coexistence line. 

While the inability to equilibrate in this domain was noted by Gibson and Wilding in 
their seminal study |30|, they did not discuss the "glass transition" (which is not a true 
transition in the thermodynamic sense). The GT is a concept useful for describing the 
manner in which viscous liquid systems fall out of equilibrium on cooling, or regain it during 
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heating. It is better described as a "glass transformation zone" within which the system 
is neither fully arrested nor fully equilibrated. It may be studied in simulation, as it is 
in the laboratory [jij, by "scanning calorimetry." In scanning calorimetry the enthalpy is 
monitored continuously as the system attempts, and increasingly succeeds, to explore all its 
degrees of freedom as temperature rises from low values where all motions except vibrations 



are frozen out 



34|. Typically the range of temperature over which the transition extends is 



the range needed to change the relaxation time by two orders of magnitude, so it depends 
on the temperature dependence of the relaxation time j^l . Being kinetic in nature, this 
transition is hysteretic, as seen in our simulations. While it is usually studied by scanning 
calorimetry, it can equally be studied by volumetric methods. 

The glass temperature T g can be defined as the point at which the uptake of configu- 
rational enthalpy begins to commence (onset T g , the value usually reported by experimen- 
talists), or the temperature at which equilibrium (ergodicity) is fully restored, T' g . Each is 
defined diagrammatically in Fig. [TQ1 the distance between the two amounting to about 25% 
of the absolute value. It is a much more diffuse phenomenon than in the laboratory where 
the width is only 5% of the absolute value (due to the increased temperature dependence of 



35|. 



the relaxation time near the laboratory T g ) 

We estimate T g and T g , as in the laboratory experiment, by plotting the derivative of the 
enthalpy (apparent specific heat), during cooling and heating the systems through GT along 
isobars slightly above the critical pressure P c , at constant cooling/heating rate q c = Qh = 
10 6 go = 10 7 K/sec [Fig. [TU]. We find that for upscans in the positively sloped coexistence line 
case (b/a = 1.72), Cp shows two well-separated peaks [Fig. ITUT a)]. The high temperature 
peak T]^, is related to the fluctuations associated with the LLPT and is used to locate 
the Widom line. The second peak (at the lower temperature T = T' g ) is an "overshoot" 
phenomenon due to ergodicity restoration kinetics. It is seen in most laboratory systems 
(but not polymers) and is not observed during cooling, (a measure of the hysteretic character 



of the k 



line) [36]. 



ass transition). The lower T g is obtained from the standard construction (dashed 



Similar well-separated TV an d T' g peaks can be observed for b/a = 1.70, and 1.68 
[Fig. [TOT b.c)]. but the temperature difference between the two peaks shrinks as the b/a 
value decreases. When b/a < 1.68, the "normal" and critical fluctuations merge because of 
the similarity in their time scales, but study of the density fluctuations as reflected in the 
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compressibility of Fig. [5] shows that indeed T c > T g , and the critical point is not suppressed 
by the kinetics of "background" enthalpy fluctuations (as the collected data in Fig. ITUl might 
imply at first sight). Rather, what is happening is that the critical fluctuations in enthalpy 
are losing their thermodynamic strength due to the vanishing of any enthalpy difference 
between the alternative phases that is dictated by the Clapeyron equation for horizontal 
coexistence lines (see Fig. [2]for b/a = 1.59). Thus, at b/a = 1.59, the apparent specific heat 
plot, notwithstanding the proximity of the critical point, is indistinguishable in character 
from that previously reported for the glass transition of the low density liquid at pressures 



well below P c in the earlier study of Xu et al. 



22]. 



Figure [TT] shows that the critical fluctuation domain becomes increasingly related to the 
slow (glassy) dynamics domain as the repulsive potential becomes steeper (second length 
scale approaches the first, more closely as b/a decreases). It is unfortunate that the increase 
of the equilibrium melting point, in the same b/a range, throws the system into conflict with 
still another, and faster kinetics, that of crystal nucleation, so that the relation between the 
first two can no longer be followed for smaller b/a. 

Just as the mixing of Lennard- Jones (LJ) particles has made possible the study of super- 
cooled and glassy states of LJ, so might the mixing of Jagla particles of different dimensions 
and attractive well depths, make possible more extended studies of the critical point/glass 
temperature relations. Note that in the glass-forming LJ mixtures there is no suggestion of 
stable domain critical points, though specific heats in excess of vibrational values do increase 
sharply as temperature decreases. 

As a final remark, it is notable that the strengths of the response functions specific heat 
and compressibility in laboratory molecular glassformers also vary in opposite directions 



as T g is approached, the case of o-terphenyl being the best documented so far [37J. The 
relationship is similar to that of the response functions maxima for the present model at 
small b/a demonstrated in Figure |3] and showing that there is a temperature interval 
near the glass transition where Cp is increasing upon cooling while the Kt decreasing, with 
the only difference that there is no stable second critical point in the laboratory case. 
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VI. SUMMARY AND CONCLUSIONS 



We have investigated the loci of the response function maxima in systems with different 
coexistence line slopes. We find that for the case of positively sloped coexistence line, the Cp 
maxima line originates from the LLCP and extends into the one-phase region as a continu- 
ation of the coexistence line, while compressibility Kt shows two maxima lines [Fig. fl2T a)]. 
One of the Kt maxima lines is related to critical fluctuations, and originates from the LLCP, 
coinciding with the Cp maxima line in the vicinity of the critical point following the Widom 
line. This offers us a method to locate the LLCP from the high temperature side by track- 
ing the Cp maxima line, instead of tracking the coexistence line from the low temperature 
side where crystallization and glass transition bring huge experimental obstacles. The other 
Kt maxima line is due to the approach to the LDL spinodal, and terminates at the LDL 
spinodal at its lowest pressure point, where the TMD line also terminates. 

As the slope of coexistence line approaches zero, Cp maxima disappears in the equilibrium 
region with T > T c [Fig. [LW b)]. However, along a constant temperature path, Cp shows 
a minimum at the critical pressure P c [Fig. This is experimentally observed in water 
for which Cp decreases with increasing pressure 19J. Hence, for a system with horizontal 
coexistence line, the LLCP can still be found by the Cp minimum as function of P at 
constant T. For Kt, both of the Kt maxima lines as functions of T are related to critical 
fluctuations, and start from the LLCP and extend symmetrically above and below P c . In 
addition, a third Kt maxima line as a function of P at constant T can be defined. All 
these three Kt maxima lines converge at the LLCP, and together with the Cp minimum 
line, can be used to locate the LLCP. Since in case of coexistence line with zero slope, the 
thermal field practically consides with temperature, this third Kt maxima line gives the 
best approximation for the Widom line. 

For negatively sloped coexistence line, the phase diagram near the critical point is similar 
to the phase diagram for the positively sloped coexistence line case reflected with respect 
to critical pressure P = P c . This follows from the linear scaling theory of the critical 
point since according to Eq. (ED K T (-<f),-0) = K T (4>,0), C P (-(f),-9) = C P (<f),6), while 
ap{— 0, — 9) = —ap((p,6), which comes from the symmetry of the functions with Xi{@)-, 
and X2{@) being even, but Xi2{9) being odd. The Cp maxima line and one of the Kt 
maxima lines both originate from the LLCP and extend into the one-phase region, while 
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overlapping with each other in the vicinity of the critical point [Fig. fl2T b)]. The second Kt 
maxima line goes above the pressure of the critical point, and terminates at the point of 
maximum pressure of the HDL spinodal, where the TMD line also terminates. 

Note that the location of the critical point with respect to the TMD line is related to the 
slope of the coexistence line. When the slope of the coexistence line is positive, the critical 
point stays outside the density anomaly region; when it is negative, the critical point is 
inside the density anomaly region; when it is horizontal, the TMD line terminates right at 
the LLCP [Fig. [12]. Indeed, when the slope is positive, the volume of the low temperature 
phase is smaller than the volume of the high temperature phase. Thus if we connect these 
two phases by the isobar with P > P c , the volume along this isobar decreases with T, the 
region above the critical point corresponds to the ap > and, because «p is a continuous 
function everywhere except at the LLCP, ap also remains positive in the one-phase region for 
pressure below the LLCP. Accordingly, the LLCP lies outside the region of density anomaly. 
Analogous considerations show that, when the slope of the coexistence line is negative, the 
LLCP remains inside the density anomaly region, as is surely the case for water. 

We note that, in all three cases, the two Kt maxima lines, merging with the Kt minima 
lines, form a loop and cross the TMD line at its highest temperature point 31]. 

By changing the parameters of the Jagla interacting potential, we can obtain systems 
with different slopes of the LLPT coexistence line. We find that, when the slope of the 
coexistence line is small, the identification of the Widom line is no longer possible by tracing 
the Cp maxima. As the slope of the coexistence line approaches zero, the Cp maxima lines 
become increasingly vertical and, when the slope of the coexistence line is horizontal, it 
cannot be observed in simulations. The study of Cp maxima is best reserved for systems 
in which the slope of the coexistence line is strongly positive or negative. However, the 
response function maxima in terms of density fluctuations are still well defined, and it is 
possible to identify the Widom line by following the loci of Kt maxima. These results are 
in good agreement with the linear scaling theory. 
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TABLE I: Renormalized parameters for modified Jagla potential 30] 
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FIG. 1: Family of modified spherically symmetric two-scale Jagla ramp potentials [Tab. H]. 
The two length scales of Jagla potential are the hard-core distance r = a and soft-core 
distance r = b. The long range cutoff is r = c. We keep the potential minimum U a con- 
stant, and its depth serves as the unit of energy Uq — \Ua\, while the hard-core poten- 
tial Ur varies. This convention is different from Ref. 30j, in which Ur is kept constant 
Ur = 0.69£7o- The discretized versions of the modified Jagla potential from b/a = 1.72 
to 1.59 are shown. We use discretization steps, n\ = 60 for a < r < b and ri2 = 20 for 
b < r < c. 
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FIG. 2: The LLCP, temperature of maximum density (TMD) line and the slope of the 
LDL-HDL coexistence line for a selection of modified Jagla potentials, (a) The critical 
point (solid symbols with error bar approximately the symbol size) and the TMD line 
(open symbols) for systems with b/a = 1.72, 1.68, 1.65, 1.62, 1.59. One can see that the 
LLCP monotonically shifts to higher pressure and lower temperature as b/a decreases. 
The density anomaly region bounded by the TMD line, expands in the T-P diagram 
with decreasing b/a, and the LLCP moves from above the TMD line towards the density 
anomaly region and for b/a = 1.59, the LLCP locates right on top of the TMD line, (b) 
The slope of the LDL-HDL coexistence lines decreases as the b/a value decreases. When 
b/a approaches 1.59, the slope decreases to zero. 
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FIG. 3: Specific heat Cp for systems with different b/a. Only equilibrated results are 
shown. The gray area indicates the T > T c region, (a-d) For b/a > 1.65, one can see that 
Cp shows maxima at pressures P > P c , and as P c is approached, the increase in Cp starts 
at lower T but becomes sharper. As the pressure increases, the Cp peak moves to higher 
T, indicating a positive slope of the Cp maxima locus, which follows the coexistence line 
for these models. For pressure P < P c , Cp monotonically increases without any maximum, 
(e) For b/a = 1.62, Cp monotonically increases without showing any peak also for P > P c , 
except at the highest pressure studied P = 0.800. The system enters a glassy state at 
lower T where no equilibrium data can be obtained, (f) For b/a = 1.59 with horizontal 
coexistence line, no Cp maxima can be found for the equilibrium states with T > T c for 
both P > P r and P < P r . 
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FIG. 4: Phase diagram with specific heat Cp for positively sloped coexistence line case 
(b/a = 1.72) and horizontal coexistence line case (b/a = 1.59). Isochores (solid gray), 
TMD line (dash purple), and LLCP (red solid circle) are shown, (a) b/a = 1.72, the lines 
of equal Cp (solid green), change from Cp = 2 far away from the LLCP to Cp = 10 close 
to the LLCP with interval ACp = 1. Cp maxima locus crosses the lines of equal Cp at 
the points of their maximal pressure, and follows the coexistence line into the one-phase 
region, then sharply turns upwards to higher pressures and becomes almost vertical, (b) 
For b/a = 1.59, the lines of equal Cp (solid green), change from Cp = 2.0 far away from 
the LLCP to Cp = 4.5 close to the LLCP with interval ACp = 0.5. No Cp maxima can 
be observed before the system either goes into glassy states or crystallizes. However, one 
notes that the Cp is symmetric with respect to the critical pressure. 
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FIG. 5: Compressibility Kt for models with different b/a. The maxima in Kt both for 
P < P c (open symbols) and for P > P c (solid symbols with arrows pointing to the peaks) 
are shown, (a-e) For b/a < 1.62,the magnitudes of Kt for P > P c are much larger than 
for P < P c . The maxima for P > P c correspond to critical fluctuations, while the maxima 
for P < P c correspond to the approach to the LDL spinodal. (f) For b/a = 1.59, with the 
horizontal coexistence line, the Kt below and above P c are almost identical, with equal 
magnitudes of their maxima. 
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FIG. 6: Phase diagram with compressibility Kt for positively sloped coexistence line case 
with b/a = 1.72 and horizontal coexistence line case with b/a = 1.59. (a) For b/a = 1.72, 
the lines of equal Kj- (solid green), changes from Kt = 6 far away from the LLCP to 
Kt = 14 close to the LLCP with interval AAV = 2. One can see that the loci of the two 
Kt maxima cross the lines of equal Kt at points of their maximal and minimal pressures, 
and are not symmetric with respect to P c . The locus with higher magnitude of Kt max- 
ima, which corresponds to the critical fluctuations, merges to the LLCP. The locus with 
the lower magnitude of Kt maxima extends below T c , and terminates at the minimum 
pressure point of the LDL spinodal (triangle symbol), where the TMD line also termi- 
nates. This brunch of Kt maxima also crosses the TMD line at the point of its maximal 
temperature, (b) For b/a = 1.59, the lines of equal Kt (solid green), changes from Kt = 2 
far away to Kt = 6 close to the LLCP with interval AKt = 1, The loci of equal Kt form 
symmetric loops around the LLCP, and both Kt maxima lines merge to the LLCP with 
equal magnitudes of their maxima. 
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FIG. 7: Cp maxima and Kt maxima obtained from the linear scaling theory of the crit- 
ical point with = 45°, 30°, 15°, 0. The coexistence line (solid black) has the slope 
dPj dT = tan 0. One can see that for more positively sloped coexistence line with larger 0, 
at P > P c , Cp maxima and Kt maxima lines have longer overlapping segments, while for 
smaller 0, the merged section is shortened, and disappears (no overlapping) when = 0, 
with a horizontal coexistence line. The second Kt maxima line at P < P c corresponding 
to the LDL spinodal can be seen for <fi = 30° and 15°. Due to the lack of metastability 
of the states in the theory, the Kt maxima at P < P c terminates at the coexistence line, 
instead in reality, it terminates at the minimum pressure point of the LDL spinodal. For 
= with horizontal coexistence line, a second Cp maxima line at P < P c merge to the 
LLCP. As we know, the linear scaling theory can predict the critical phenomena only close 
to the LLCP, a closer study near the LLCP will follow. Further, the TMD line (dash) is 
also shown, which is the locus of otp = 0. One notices that as the slope of coexistence line 
goes from positive to horizontal, the TMD line shifts from being below the LLCP, to right 
on top of it, which is consistent with our results in DMD simulation. Here a = 0.47, and 
k = 0.47. 
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FIG. 8: Cp maxima and maxima loci on T-P plane obtained from the linear scaling 
theory of the critical point, with small <fi = 8°, 4°, 6°, 0. One can see that as we decrease <fr, 
with the slope of coexistence line approaching zero, the Kt maxima locus changes are not 
significant, while the Cp maxima become more and more vertical. For a horizontal slope, 
two Cp maxima lines converge at T < T c and enter the LLCP horizontally from below. 
There is no overlapping between Kt maxima and Cp maxima lines. Here we use a = 2, 
and k — 1, which is good to present the small or zero coexistence line case. 
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FIG. 9: Comparison of Cp and Kt from linear scaling theory of the critical point. (a,b) 
Cp and Kt for the positively sloped coexistence line with = 30°, respectively. The lines 
of equal Cp and Kt are shown. Both Cp and Kt contours loop around LLCP, with the 
locus of the maxima follow the coexistence line into the one-phase region. (c,d) Cp and 
Kt for the horizontal coexistence line with = 0, respectively. One can see that Cp and 
Kt are both symmetric below and above P c . For clarity, (c) displays a narrow vicinity of 
the LLCP. Here we choose a = 0.47, and k = 0.47. 
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FIG. 10: The comparison of Cp upon slow cooling and heating for a selection of modi- 
fied Jagla models. The cooling/heating rate is q c = Qh = 10~ 6 <7o- An isobaric path is 
selected at P closely above P c in the one-phase region. The equilibrium values of Cp are 
also plotted for comparison. The gray area is the temperature region for T > T c . (a- 
c) For b/a = 1.72, 1.70 and 1.68, two peaks of Cp upon heating can be observed. The 
high temperature peak TV, arising from the presence of the Widom line, is related to the 
LLPT. The low temperature peak upon heating, T g , corresponds to the ergodicity restora- 
tion slightly above the glass transition T g obtained from standard construction (dash line). 
One can see that the distance between the glass transition peak T' g /T g and TV decreases 
as the b/a value decreases, with the LLCP moving closer to the glass transition, (d-f) For 
b/a = 1.65, 1.62 and 1.59, upon heating, only one peak of Cp can be observed, and this 
peak shifts below T c , while for other cases it is well above T c . One also notices that system 
shows a prominent peak in the studied P region (Fig. 5). In these models with small 
coexistence line slope, the enthalpy fluctuations play less role in the critical fluctuations, 
while the density fluctuations is the leading term. And as the LLCP being pushed closer 
to T g , the critical fluctuations become suppressed by the glass transition, where Cp only 
shows one peak upon heating. 
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FIG. 11: Relative positions of the glass transition temperature T g (onset), T' g (upper 
limit), critical temperature T c , and the locus of Cp and Kt maxima near the LLCP, for 
models of different b/a. Only the high pressure branch of the Kt maxima locus is shown. 
The LLCP shifts to lower temperature and higher pressure as b/a decreases. Tg follows 
the same trend of T c , but the temperature difference (gray area) between Tg and T c de- 
creases as b/a decreases. The glass transformation range (between T g and T' g ) separat- 
ing glass from liquid is shown (hashed area). For b/a < 1.65, T c gets pushed close to 
Tg, where in isobaric cooling/heating, only one peak of Cp can be found, instead of two 
well-separated maxima (b/a = 1.72 — 1.68). The locus of Cp maxima increasingly sep- 
arates from the locus of the Kt maxima with decreasing b/a, and is no longer seen for 
b/a < 1.65. The Kt maxima line is a better representative of the Widome line for the case 
of small and zero sloped coexistence line. T c lies in the ergodic domain for the b/a stud- 
ied. For b/a < 1.59, the LLCP cannot be obtained in our study by equilibrium molecular 
dynamics. 
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FIG. 12: The relationship between the slope of the coexistence line and the behavior of 
the Widom line and the TMD line, (a) For the case of positively sloped coexistence line, 
the Cp maxima line originates from the LLCP and extends into the one-phase region as 
a continuation of the coexistence line, while the Kt shows two maxima lines. One Kt 
maxima line is related to critical fluctuations, and originates from the LLCP, coinciding 
with the Cp maxima line in the vicinity of the critical point forming the Widom line. The 
other Kt maxima line corresponds to the approach to the LDL spinodal, which terminates 
at the LDL spinodal at its lowest pressure point, where the TMD line also terminates (tri- 
angle symbol), (b) For the case of negatively sloped coexistence line, the phase diagram 
is similar to the mirror image of the positively sloped coexistence line case, with Cp and 
Kt maxima lines originating from the LLCP and extending into the one-phase region, 
and overlapping in the vicinity of the critical point. The second Kt maxima line goes 
above the pressure of the critical point, and terminates at the HDL spinodal, where the 
TMD line also terminates. We note that in both cases, both Kt extrema lines form a loop 
in the T-P plane, which intersects the TMD line at its maximum temperature point of 
the TMD line (not shown). In case (a) the critical point lies outside the density anomaly 
region while in case (b) it lies inside the density anomaly region. Both graphs are con- 
structed using the linear scaling theory of critical point with a — k — 0.47, and <ft — 25° 
and —25°, respectively. The lines of the spinodals are drawn as interpolation between the 
critical point and the extrapolated crossing point of the TMD and Kt maxima lines be- 
yond the coexistence line. 
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FIG. 13: The behavior of Cp and Kt along different isotherms as functions of P in sim- 
ulations with a horizontal coexistence line b/a = 1.59 (a,c), and the linear scaling the- 
ory with = (b,d). (a) Cp in simulations, T changes from 0.276 to 0.442 with interval 
AT = 0.180, above T c = 0.226. (b) Cp from linear scaling theory. Both results agree well 
with each other, that Cp shows a minimum at P c for all temperatures, and as T — > T c , the 
minimum value of Cp increases, and the valley of the minimum gets narrower. This offers 
a way to track the critical point by isotherms at equilibrium temperatures, instead of iso- 
bar cooling into lower temperature, (c) Kt in simulations, T changes from 0.239 to 0.369 
with interval AT = 0.180, above T c = 0.226. (d) Kt from linear scaling theory. Again, the 
theory is consistent with the simulation results, that showing Kt has maxima at P c for 
all temperatures and, as T — > T c , the maximum value of Kt increases and the peak gets 
sharper. 
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